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Abstract 

We study dynamics of populations of resonantly coupled oscillators having dif- 
ferent frequencies. Starting from the coupled van der Pol equations we derive 
the Kuramoto-type phase model for the situation, where the natural frequencies 
of two interacting subpopulations are in relation 2:1. Depending on the param- 
eter of coupling, ensembles can demonstrate fully synchronous clusters, partial 
synchrony (only one subpopulation synchronizes), or asynchrony in both sub- 
populations. Theoretical description of the dynamics based on the Watanabe- 
Strogatz approach is developed. 
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1. Introduction 

Models of coupled autonomous oscillators are used to describe synchro- 
nizationphenomena [l|, 0] appearing in many physical 0, 0, [H] and biological 
d, 0, H systems. In the case of a weak coupling, a phase model description 
is appropriatejeading to the famous Kuramoto model 13, HI and its modi- 



fications [12], [l^l . One of the main assumption behind the derivation of these 
models is that the oscillators are in resonance, i.e. their frequencies are close 
to each other (even when a bimodal distribution of frequencies is considered 



(see [IJ, [15| and references there) , one assumes that the distance between the 
peaks is small). Recently, we considered ensembles of oscillators consisting of 
non-resonantly coupled groups [lij . i.e. those with frequencies that are far from 
each other and far from resonances. 

In this paper we study synchronization effects in ensembles where different 
groups of oscillators are in a non-trivial resonance relation 2:1. First, we 
derive general equations describing these interacting subpopulations in the phase 
approximation. Then we demonstrate numerically regimes of complete and 
partial synchrony (in the latter case one subpopulation synchronized while the 
other not). Furthermore, we develop a theory based on the Watanabe-Strogatz 



approach [17|, |18|, |19|, |20|, |2lJ that allows one a description in terms of dynamical 



equations for the order parameters. 
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2. Basic Model 



While typically one considers ensembles of oscillators with close frequencies 



that interact resonantly here we focus on the two populations having nat- 
ural frequencies lo and = 2lo. Therefore, we spend more space than usual 
describing the derivation of the phase model. We start with two coupled van 
der Pol oscillators 

X - ^J.i{l - x^)x + oj'^x = fi(x,x,y,y) , 

y - M2(i - y'^)y + ^'^y = f2{x, x, y, y) , 

and following a standard procedure write the averaged equations for the slow 
varying complex amplitudes A{t) = e~'"*(a; — ixjiS)^ B{t) = e~'^^*(a; — ix/U): 

A=^{1- \A\y4)A--{h{x,x,y,y)e-^-')) , 

B^^il- \B\y4)B - l(/2(x,i,y,y)e-"*)) . 

The interacting terms that survive the averaging are those with dependence 
/i ~ e'"* and /2 ~ e'"*. Thus, due to the resonance condition 2uj — 57, /i 
should contain a product of x and y, while /2 should contain the square of 
x. Therefore the simplest polynomial terms that yield a coupling between two 
oscillators are /i = cixy + C2xy + c^xy + c^xy and /2 = dix^ + d2xx + d^x^ . 
Correspondingly, the averaged equations can be written as 

i = ^(1 - \A\^/A)A + aie'"'A*B , 
B = ^il-\B\y4)B + a2e''''A\ 

with some complex coupling constants cri,2e*"^'^ that can be expressed in terms 
of constants ci_4,c?i_3. As the next step, we use smallness of cri.2 compared to 
/xi.2, so that the deviations of the amplitudes \A\, \B\ from the limit cycle values 
\A\ ^ \B\ =2 are small. Then, substituting A = 2e*'^(*) and B = 2e*''^(*), we 
obtain for the phase dynamics 

(j) = 2a\ sin(?/; — 20 + a\) , 
-0 = 2a-2 sin(20 - ?A + 02) . 

By shifting one of the phases '\\} — '\\)' — a\ and rescaling the time variable 
2(tTi + a-2)t = t' we can reduce the dynamics to a system with two parameters 
[1 = a\l (cTi +(J-i) and 7 = ai+a2 only (we use the same letters for new variables) 

dcj) . , . 
— = ^ism(tp ~ 20) , 

^ = (l-M)sin(20-V + 7)- 
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These equations describe two coupled phase oscillators cj) and ip with a relative 
coupling strength ^ G [0, 1]. Because of symmetry 0, 7 ^ —ip, —tp, —7 we vary 
the phase shift in the range 7 G [0, tt]. The parameters /i and 7 are functions of 
the coupling terms in Eqs. (P). 

Now we generalize to ensembles of oscillators, assuming that each unit in 
a population with single frequency oj interacts with every unit in a population 
with double frequency J7 = 2a;: 

^^^}:^sm(20,-^.+7). 

here N^^, Nq are the sizes of subpopulations. Equations ([2]) describe the basic 
model that we will investigate in the following. It consists of two groups of 
oscillators with a frequency ratio 2:1. Each group is composed of identical 
oscillators. One oscillator of a group is coupled to all oscillators of the other 
group, and vice verse. We assume that there is no interaction within one group. 



3. Dynamical regimes and their characterization 

We first present numerical results of simulations of ensemble ([2]) , setting the 
relative coupling strength /i = 0.5 and the number of oscillators to be equal in 
each group Nfi = N^^ = N . Our main attention here is to the dependence of 
the dynamics on the phase shift 7 and on different initial conditions. 

We illustrate a nontrivial regime of the interaction of two populations in 
Fig. [TJ Here, for N — 200 and 7 = 2.8, by integrating Eqs. ([2]) we observe that 
single- frequency oscillators {(f>k, blue dashed curves) form two clusters that differ 
by TT, while double- frequency oscillators {ipk, red full curves) remain distributed 
in some range of phases. 

To characterize the synchronization properties and the clustering, we adopt 



the Daido generalized order parameters [12, l22|, l23| , calculated separately for 
double- and single-frequency ensembles: 

^1 ^ ]^ E [^j'-^fe (^)] ' (^) = 1^ E e'^p t^j'^'' (*)] ■ (3) 
fe=i fc=i 

The physical meaning of these quantities is clear from considering the case of 
large ensembles, then it follows from ([3]) that Zj, Yj are the j-th Fourier modes 
of the distributions of the phases. While the usual Kuramoto order parameters 
Z\ , Y\ are suitable for characterization of distributions having a single maximum 
(single clusters), the second order parameters Z2, Yi allow us to reveal 2-cluster 
states (distributions with two humps with the phase difference tt) - at these 
states these parameters have absolute value one, while the Kuramoto parameters 
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(a) time (b) time 



Figure 1; (a) Phases of coupled ensembles l(2} as functions of time for 7 = 2.8, = 200, and 
/i = 0.5 and the corresponding evolution of the generalized order parameter (b). In (a) red full 
lines are rf>k (only 10 phases from the population are shown for clarity) and blue dashed lines 
are tp^.. In (b) red full lines, blue dashed lines and green dotted lines correspond to j = 1, 2, 3. 

Zi , Yi vanish. Similarly, Z3 , Y3 are suitable for revealing an existence of three 
clusters with phase shift ^ (or, more generally, three- hump distributions). 

In Fig. (HJb) the evolution of the order parameters for the dynamics depicted 
in Fig. [TJa) is shown. After a transient time interval, one can see that IZ2I = 1 
while \Zi \ < 1 and jZaj < 1, as expected. We note here that in the case IZ2I = 1, 
the exact values of order parameters {Zi^^l are irrelevant as they characterize 
the partition between two clusters seen in Fig. [Ija). As these clusters differ 
by TT, this partition has no effect on the dynamics, where only the values 2(f>k 
entry (cf. Eq. The variations of the order parameters lYj-j, j = 1,2,3 

characterize the oscillating distribution of the phases of the double-frequency 
oscillators. 

Next, we want to characterize the dynamics of the order parameters in de- 
pendence on the parameter 7. Therefore, after a certain transient time we calcu- 
lated the minimal and maximal values of the amplitudes of the order parameters 
in course of their evolution. Hence, a finite interval between the maximum and 
the minimum characterizes a range of variations of the order parameter: e.g., 
for the regime presented in Fig. [TJ variations of \Yj\ are finite while there is no 
variations in \Zj\. 

We have found that dynamical regimes strongly depend on initial conditions. 
In Figs. l2l3l we show the variations of the order parameters for two sets of initial 
conditions: in Fig.[2]both phases are uniformly distributed in the interval [0, it); 
in another set Fig. [3] the initial conditions are chosen specially, according to the 
theory we develop below in Section H) As discussed above, the most relevant 
are the data for Yi and Z2, while values of Zi,Z^ giving a partition between 
subclusters of single-phase oscillators are irrelevant for the dynamics. For < 
7 < ~ 2.094 both these order parameters are one, what means that in both 
subpopulations full synchrony establishes. For 7ci < 7 < tt the order parameter 
|Yi| < 1 what means asynchrony in the double- frequency subpopulation. For 
initial conditions in Fig. |3] the single-frequency subpopulation is synchronized 
in this range, while for the setup of Fig. [5] the single-frequency population is 
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Figure 2: Variations of the order parameter \Zj{-y)\, |5^(7)| for Nq = Noj = 200, /i = 0.5 in 
dependence on 7. Here the initial distribution is a uniform one over half of the phase circle: 

< (pk,i'k < T- 




Figure 3: The same as in fig. ((2]|, but with a different initial condition; The oscillators were 
initially transformed according to the WS theory (see Eq.lJjl below). 

synchronous up to jc2 ~ 2.98 and asynchronous for 7 > jc2- Thus, the coupled 
ensembles demonstrate regimes of full synchrony for 7 < 7ci , partial synchrony 
for 7ci < 7 < 7c2 and asynchrony (unless special initial conditions are chosen) 
for 7c2 < 7 < TT. 

To reveal the type of dynamics in regimes of partial synchrony and asyn- 
chrony, we show in Fig. |4] the two-dimensional projections of trajectories on the 
planes of main order parameters. One can see that in both cases the dynamics 
is two- frequency quasiperiodic; this is also confirmed by calculations of Poincare 
maps, on which the attractor form one-dimensional lines. 



4. Watanabe-Strogatz ansatz 

In this section we apply the WS ansatz allowing us to describe the system 
of coupled phase oscillators with a few global variables. In this way one is able 
to analyze such systems analytically. We begin with a sketch of the WS theory 
according to 2l|, for an original formulation see 17, 18]. One starts with an 
ensemble of identical phase oscillators with frequencies Ld{t) driven by a force 
H{t) according to 

= uj{t) + Im[H{t)e-'^''] , k = 1, iV > 3 , 
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Figure 4: Projections of the dynamical regimes on the planes of order parameters, for A'^ = 200 
and /i = 0.5. (a): 7 = 2.8, here IZ2I = 1 (blue circle) while Yi varies in some range (red curve). 
(b),(c): For 7 = 3.08 both order parameters vary quasiperiodically. 



and performs a transformation to new microscopic variables i3fe and gfobal vari- 
ables z,( (z is complex, is real) according to 

(4) 



with additional conditions ^ 6*'''= — and Re(^ e*^''*) — 0. Then, the new mi- 
croscopic phases 1?^ are constants of motions provided the macroscopic variables 
z, ( satisfy the WS equations 

(77 1 

-^^u;{t) + -{H{t)-z'H*{t)) , 

^=u;it)+lm{z*Hit)). 

As discussed in pH . the complex variable z is roughly proportional (but not 
exactly equal) to the Kuramoto order parameter of the population, while C is a 
shift between individual phases and the phase of z. Indeed, from the transfor- 
mation (U]) it follows that 

(^'^)-]^E i + ,.e-(..K) -<^'^)(-'C,^.). (5) 

Only in the case of a uniform distribution of constants on the interval [0, 2tt) 
and in the thermodynamic limit TV — >■ cxi, the variable z coincides with the Ku- 
ramoto order parameter, because the dependence on ( and "dk hi ([5]) disappears 
and z = {e"^). 

To apply the WS theory to our system 1^ we perform a transformation 
20fc = 0k and rewrite the system as 



V^fe - (1 - /i)(sin {0j ~i^k+ 7)) - (1 - M)Im(Z2e-*'^'=e^^) 



(6) 



Now we are able to apply the WS transform with two sets of WS variables zi^2, 
C1.2, for subpopulations of oscillations having single and double frequencies. 
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correspondingly. As discussed above, the consideration extremely simplifies in 
the case of a uniform distribution of the corresponding constants of motion "dk 
and in the thermodynamic limit, where a simple representation of the order 
parameters via the WS variables holds: 

zi = Z2 and Z2 = Yi . (7) 

In this case, the closed WS system reads 

at /g\ 

C!£2 _ (1^ 2 

dt ~ 2 ^'^^ '^'^^ > ■ 
With introduction of amplitudes pi_2 = |zi,2| and the phase difference 4* = 
arg(z2) — arg(zi) we obtain a three-dimensional system 

Pi = (1 - pI) P2 cos {^) , 

P2 = ^ (1 - A*) Pi (1 - pD cos - 7) , (9) 

^ = pi ^-^sinfvlf - 7) - ^ ^p2sin\I' . 

2 P2 Pi 

We first discuss the steady states in system © and their stability. The 
steady state corresponding to a full synchrony is 

p(°).p(«)=l, tanvl;(°)^ ^ (M-l)sin7 _ ^ 
'^^ ' (/^- l)cos7-2^ ' 

It is stable if cos 7 > max(i^j^, ■^). (There is another state with = vl/(o) + vr 
which is unstable). 

Two stable steady states with partial synchrony are possible. For parameter 
values —1 < C0S7 < p > 1/3 the state 



V 4^ cos 7 + 1 — /i 



is stable while for — 1 < cos 7 < /i < 1/3 another state 



V [p. -I) COS 7 - /i 

is stable. The fully asynchronous state pi = p2 = is unstable for 7 < tt. 

The special case 7 = tt deserves a separate analysis. In this case the dynamics 
is described by equations 

pi =m(1-P?)p2COS* , (11) 
P2 = i(l-p)pi(l-p^)cos*, (12) 

^=(i^Pii±^-M^P2)smvI,. (13) 
V 2 p2 Pi J 
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One can see that this system has two integrals. One is obtained by dividing 
(fTTj) and p2|) and integration; another one is obtained by first using the first 
integral to express P2(pi), and then dividing (ITT]) and whith subsequent 
integration. Thus, the resulting dynamics is conservative and periodic. 

The analytical results above explain Fig. [31 where the initial conditions have 
been chosen according to Eq. Q with a uniform distribution of the constants z?fc . 
In Fig. |3] one observes a regime of full synchrony corresponding to steady state 
^(0)^^(0) foj. ^ < = 27r/3 and the steady state p^^^,^'^^) for 7 > 27r/3, in 
according with relations above and /i = 1/2. For 7 = tt one observes oscillations 
of the order parameters. 

In order to explain Fig. [2] we have to go beyond the assumption used at 
the derivation of ([5]), namely of the uniform distribution of constants dk and 
of thermodynamic limit. In general case, instead of (O we have to use (O for 
Z2 — ^2(^1, Cii '^i^') and Yi = Yi{z2, C2, ^i^^)- Now the full system of equations 



The fully synchronous state where both populations are completely synchronized 
is the same as in system because as one can see from ([S]), for \z\ = 1 we 
again have {e^^) = z, but all other states are generally different. In particular, 
partially synchronous regimes are not steady states but quasiperiodic ones as in 
Figs. I1I2|4I We discuss their stability in the next section. 

5. Stability of partially synchronous state 

The WS theory above allows us to describe analytically the transition full 
synchrony — > partial synchrony as loss of stability of the fully synchronous 
state ([in]), but the analysis of the transition from partial synchrony (where one 
ensemble forms a synchronous cluster while other one is asynchronous) is more 
involved as it deals with quasiperiodic states. Therefore we apply here a direct 
numerical method for determining stability of clusters - calculation of so-called 
evaporation Lyapunov exponents [3, H^. We assume, according to numerics, 
that the single- frequency oscillators are different, while the double- frequency 
oscillators form a synchronous cluster, i.e. in Eq. ^ ipi — ip2 — ■ ■ ■ — i^Na — ^• 
Then, in the limit iVo — > 00, the deviation of one element from the cluster is 
governed by 



IS 




M {Yi - zlY*) 
2/iIm(^i*Yi(i)) 



(14) 



(1 - ii)ln\{zlZ2) 




8 








3.08 



3.12 



-1.5 



2.2 



2.4 



2.6 



2.8 



7 



Figure 5: Evaporation exponent calculated for Af^j = 200, fi = 0.5, in dependence on 7. Initial 
conditions: phases ij>k uniformly distributed in [0,7r). The inset shows region around 7 = tt. 
For 7 < 3.05 the exponent is large and negative, while for 7 > 3.1 the exponent vanishes. 



Thus, the growth rate of 6ip is determined by the evaporation exponent 



The resuhs of calculation of this exponent are presented in Fig. [Sj One can 
see that the synchronized cluster of double-frequency oscillators is strongly sta- 
ble for 7 < 3.1 while for 7 > 3.1 the evaporation exponent vanishes. This 
means marginal stability of the synchronized state, what is consistent with the 
observation of non-synchronous dynamics for appropriate initial conditions. 

6. Conclusion 

In this paper we have studied a novel model of resonantly interacting multi- 
frequency oscillator populations. As the simplest setup we have chosen a situ- 
ation where oscillators are divided in two subpopulations; some have natural 
frequency w while other ones have the double frequency D, — 2lu. Such a setup 
can be easily generalized to a general case of two subpopulations in a reso- 
nance D, : Lo = m : n. Moreover, one could study many subpopulations having 
resonantly related frequencies uji : ll>2 '■ ^3 '■ ■ ■ ■ = mi : TO2 : 7713 : • • • , the 
generalization of equations ([2]) to this case is straightforward. 

Our main finding is that depending on the parameters of the coupling of 
two ensembles, one observes regimes of full synchrony (both subpopulations 
form fully synchronous clusters), partial synchrony (one subpopulation syn- 
chronized while other is asynchronous) and no synchrony (both subpopulations 
asynchronous). The latter two regimes demonstrate quasiperiodic dynamics. 
To analyse these regimes we applied the Watanabe-Strogatz theory and derived 
the equations for macroscopic variables describing distributions of oscillators in 
subpopulations. This allowed us to identify the transition from full to partial 
synchrony as a transcritical bifurcation of this system. To characterize the tran- 
sition from partial synchrony to asynchrony we used the method of evaporation 
Lyapunov exponents. 




d 1-iJ 



sin(20j -^p + -/) 
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In this paper we restricted our attention to the case of identical oscilla- 
tors in subpopulations. For the study of non-identical ensembles the powerful 
Ott-Antonsen theory [2^ 27 1 can be adopted, these results will be reported 



elsewhere. 
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